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Recent experimental progress lias made it possible to detect in real-time single electrons tunnel- 
ing through Coulomb blockade nanostructures, thereby allowing for precise measurements of the 
statistical distribution of the number of transferred charges, the so-called full counting statistics. 
These experimental advances call for a solid theoretical platform for equally accurate calculations 
of distribution functions and their cumulants. Here we develop a general framework for calcu- 
lating zero-frequency current cumulants of arbitrary orders for transport through nanostructures 
with strong Coulomb interactions. Our recursive method can treat systems with many states as 
well as non-Markovian dynamics. We illustrate our approach with three examples of current ex- 
perimental relevance: bunching transport through a two-level quantum dot, transport through a 
nano-electromechanical system with dynamical Franck-Condon blockade, and transport through co- 
herently coupled quantum dots embedded in a dissipative environment. We discuss properties of 
high-order cumulants as well as possible subtleties associated with non-Markovian dynamics. 

PACS numbers: 02.50.Ey, 03.65.Yz, 72.70.+m, 73.23.Hk 



I. INTRODUCTION 

Electron transport through nanoscale structures is a 
stochastic process due to the randomness of the individ- 
ual tunneling events. Quantum correlations and electron- 
electron interactions can strongly influence the transport 
process and thus the statistics of transferred charges. 
Full counting statisticsEHj concerns the distribution of 
the number of transferred charge, or equivalently, all cor- 
responding cumulants (irreducible moments) of the dis- 
tribution. Conventional transport measurements have 
focused on the first cumulant, the mean current, au' 
in some cases also the second cumulant, the noise. 
Higher order cumulants, however, reveal additional infor- 
mation concerning a variety of physical phenomena, in- 
cluding quantum coherence, entanglement, disorder, and 
dissipationjj For example, non-zero higher-order cumu- 
lants reflect non-Gaussian behavior. Counting statistics 
in mesoscopic physics has been a subject of intensive the- 
oretical interest for almost two decades, but recently it 
has also gained considefaiJ|e experimental interest: in 
a series of experiments,ErtII high order cumulants and 
even the entire distribution function of transferred charge 
have been measured, clearly demonstrating that count- 
ing statistics now has become an important concept also 
in experimental physics. 

The theory of counting statistics was first formu- 
lated by Levitov and Lesovik for jMn-interacting elec- 
trons using a scattering formalism.Eja Subsequent works 
have focused r-G|p--| the inclusion of interaction effects 
in the theoryJlallj In one approach. Coulomb interac- 



tions are incorporated via Markovian (generalized) mas- 
ter equatwpns as originally developed by Bagrets and 
Nazarov.c2l This approach is often convenient when con- 
sidering systems with strong interactions, e.g.. Coulomb- 
blockade structures. More recent developments-include 
theories for finite- frequeiicv counting statistics j£d condi- 
tional counting statistics,E3 connections ±D_entanglement 
entropyo and to fluctuation theoremsol£3 aad jfixten- 
sions to systems with non-Markovian dynamics. EZtl^ The 
last topic forms the central theme of this paper. 

We have recentU^ jMiblished a series of papers on 
counting statistics. E3c3 Previous methods for evaluat- 
ing the counting statistics of systems described by mas- 
ter equations had in practice been limited to systems 
with only a few states, and in Ref. ^ we thus devel- 
oped techniques for calculating the first few cumulants 
of Markovian systems with map¥ states, for example 
nano-electromechanical systems.E^ In Ref. |2^, Braggio 
and co-workers generalized the approach by Bagrets and 
Nazarov by including non-Markovian effects that may 
arise for example when the coupling to the electronic 
leads is not weak. The methods presented in these papers 
were subsequently unified and extended in Ref. where 
we presented a general approach to calculations of cumu- 
lants of arbitrary order for systems with many states as 
well as with non-Markovian dynamics. The aim of the 
present paper is to provide a detailed derivation and de- 
scription of this methofL-which recently has been used 
in a number of worksJ^J^Zl as well as to illustrate its use 
with three examples of current experimental relevance. 

The paper is organized as follows: In Sec. || we in- 
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troduce the generic non-Markovian generahzed master 
equation (GME) which is the starting point of this work. 
The GME describes the evolution of the reduced den- 
sity matrix of the system, which has been resolved with 
respect to the number of transferred particles. Memory 
effects due to the coupling to the environment as well as 
initial system-environment correlations are included in 
the GME. Within this framework it is possible to calcu- 
late thfiJinite-frequency current noise for non-Markovian 
GMEsEEI as we will discuss in future works. Section || 
concludes with details of the superoperator notation used 
throughout the paper. 



In Sec. Ill we develop a theory for the zero-frequency 
cumulants of the current. The cumulant generating func- 
tion (CGF) is determined by a single dominating pole 
of the resolvent of the memory kernel, and its deriva- 
tives with respect to the counting field evaluated at zero 
counting field yield the cumulants of the current. Even 
in the Markovian case it is difficult to determine analyt- 
ically the dominating pole and in many cases one would 
have to find it numerically. Numerical differentiation, 
however, is notoriously unstable, and often one can only 
obtain accurate results for the first few derivatives with 
respect to the counting field, i. e. the cumulants. In order 
to circumvent this problem, we develop a numerically sta- 
ble recursive scheme based on a perturbation expansion 
in the counting field. The scheme enables calculations 
of zero- frequency current cumulants of very high orders, 
also for non-Markovian systems. Some notes on the eval- 
uation of the cumulants are presented, with the more 
technical numerical details deferred to App. |^ 

Section gives a discussion of the generic behavior 
of higlt£)rdcr cumulants. As some of us have recently 
shown,0 the high-order cumulants for basically any sys- 
tem (with or without memory effects) are expected to 
grow factorially in magnitude with the cumulant order 
and oscillate as functions of essentially any parameter as 
well as of the cumulant order. We describe the theory 
behind this prediction which is subsequently illustrated 
with examples in Sec. 

Section ^ is devoted to two Markovian transport mod- 
els of current research interest, which we use to illustrate 
our recursive scheme and the generic behavior of high- 
order cumulants discussed in Sec. [V. We start with a 
model of transport- 



veloped by Belzig.E 



Jirough a two-level quantum dot de- 
3 Due to the relatively simple ana- 
lytic structure of the model, it is possible to write down 
a closed-form expression for the CGF, allowing us to de- 
velop a thorough understanding of the behavior of high- 
order cumulants obtained using our recursive scheme. We 
study the large deviation function of the system,L3 which 
describes the tails of the distribution of measurable cur- 
rents, and discuss how it is related to the cumulants. 

The second example concerns charge transport cou- 
pled to quantized mechanical vibrations as considered 
in a recent series-, of papers on transport through 
single rrp]ecjLilcsc2rE3 and other nano-electromechanical 



ticipating in transport the matrix representations of the 
involved operators are of large dimensions and it is nec- 
essary to resort to numerics. We demonstrate the numer- 
ical stability of our recursive algorithm up to very high 
cumulant orders (~ 100) and show how oscillations of the 
cumulants can be used to extract information about the 
analytic structure of the cumulant generating function. 
We calculate the large deviation function and show that 
it is highly sensitive to the damping of the vibrational 
mode. 



systems! 



)ue to the many oscillator states par- 



Section VI concerns the counting statistics of non- 
Markovian systems. We consider a model of non- 
Markovian electron transport through a Coulomb- 
blockade double quantum dot embedded in a dissipative 
heat bath and coupled to electronic leads. The dynam- 
ics of the charge populations of the double dot can be 
described using a non-Markovian GME whose detailed 
derivation is presented in App. |c[ We study the behav- 
ior of the first three cumulants thus extending jxrevious 
studies that have been restricted to the noise.HaWe fo- 
cus in particular on the infiuence of decoherenceE3 on the 
charge transport statistics. Finally, we discuss possible 
subtleties associated with non-Markovian dynamics and 
we provide the reader with a unifying point of view on a 
number of results reported in previous studies as well as 
in the examples discussed in this paper. 

Our conclusions are stated in Sec. VU. Appendix ^ 
describes the numerical algorithms used to solve the re- 
cursive equations for high-order cumulants, while Apps. 
P and |c| give detailed derivations of the Markovian GME 
for the vibrating molecule and the non-Markovian GME 
for the double-dot system, respectively. 



II. GENERALIZED MASTER EQUATION 

The generic transport setup under consideration in this 
work is depicted in Fig. ^ A nanoscopic quantum sys- 
tem is connected by tunneling barriers to two electronic 
leads, allowing for charge and energy exchange with the 
leads. Typically, the quantum system consists of a dis- 
crete set of (many-body) quantum states. Moreover, the 
system is coupled to an external heat bath to and from 
which energy can fiow. We will be considering a trans- 
port configuration, where a bias difference between the 
leads drives electrons through the system. 

The quantum system is completely described by its 
(reduced) density matrix p(t), obtained by tracing out 
the environmental degrees of freedoms, i. e. the electronic 
leads and the heat bath. It is, however, advantageous to 
resolve p{t) into the components p{n, i), corresponding to 
the number of electrons n that have tumiclcd through the 
system during the time span [0,t].L£rl£j The n-resolved 
density matrix allows us to study the statistics of the 
number of transferred charges, similarly to well-known 
techniques from quantum opticsEjO We note that the 
(un-resolved) density matrix can always be recovered by 
summing over n, p(t) = X]nP('^iO- For bi-directional 
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FIG. 1: Generic transport setup. A quantum system is con- 
nected to electronic leads and a heat bath. A bias difference 
between the leads drives electrons through the system, which 
can exchange energy with the surrounding heat bath. The 
system is described by the n-resolved density matrix p(n, t) 
(see text), where n is the number of electrons that have been 
collected in the right lead during the time span [0,t]. The 
probability distribution of n is denoted as P{n,t). 



processes, the number of tunneled electrons n can be both 
positive and negative. 

The major focus in the literaluxe has been on systems 
obeying Markovian dynamics ;I3xj however recent years 
have witnessed an. increased interest in non-Markovian 
processes as well.lI3c3 In this spirit we consider a generic 
non-Markovian generalized master equation (GME) of 
the form 



dt'W{n -n',t- t')p{n', t') + 7(71, i), 

(1) 

obtained by tracing out the electronic leads and the heat 
bath. An equation of this type arises for example in the. 
partitioning scheme devised by Nakajima and Zwanzig,Ej 
and in the real-time diagrammatic technique for the dy- 
namics of the reduced density matrix on the Keldysh 
contour as described in Refs. p7p^ , |8^ -|8^. The memory 
kernel W accounts for the dynamics of the system taking 
into account the influence of the degrees of freedom that 
have been projected out, e. g. the electronic leads and 
the heat bath. Here we have assumed that the system is 
not explicitly driven by any time- varying fields, such that 
the kernel W only depends on the time difference t — t' . 
Additionally, we assume that the number of electrons n 
that have been collected in the right lead does not af- 
fect the system dynamics, and the kernel consequently 
only depends on the difference n — n' . The generic non- 
Markovian GME also contains the inhomogeneity 7(71, t) 
which accounts foe jiaitial correlations between system 
and environment. ljH Typically, W{n,t) and j{n,t) de- 
cay on comparable time scales, and 7(71, t) thus vanishes 
in the long-time limit of Eq. (Q) . 

At this point, we note that while our method for ex- 
tracting cumulants works for any GME which satisfies 
certain, rather general conditions specified in detail in 



Sec. Ill, the physical meaningfulness of the results never- 
theless depends crucially on a consistent der ivatio n of the 
n-resolved memory kernel yV(n,t). In Sec. VI B we dis- 



cuss various subtleties associated with a proper deriva- 
tion of the memory kernel for non-Markovian systems. 
In this work we use the notion of the Markovian limit 
of a general non-Markovian GME in a somewhat loose 
manner, namely by referring to the "Markovian" limit 
of Eq. (§) as the case, where W{n,t) = 'W{n)6{t) and 
j{n,t) = 0. We use this terminology for the ease of 
notation, although we are aware that the proper Marko- 
vian limit under certain circumstances may actually be 
different. For an example of this, we refer the reader to 
Ref . where it is demonstrated that the correct Marko- 
vian limit for weak coupling theories should be performed 
in the interaction picture. Since this procedure only in- 
fluences the off-diagonal elements in the weak coupling 
regime, we ignore this subtlety in the rest of the paper 
as we will not be considering such cases. In relevant 
situations this difference should be taken into account 
— it would, however, only lead to a reinterpretation of 
the non-Markovian corrections studied in Sec. VI B. The 
issue of non-Markovian behavior, its nature and distinc- 
tion from M ar i kovi an approximations, is a nontrivial and 
timely topicBEaEj which we only touch upon briefly in 
this work, but our formalism paves the way for system- 
atic studies of such problems in the context of electronic 
noise and counting statistics, for example as in Ref. 5ll 



A. Counting statistics 

In the following we introduce the notion of cumulants 
of the charge transfer probability distribution, and derive 
a formal expression for the cumulant generating function 
CGF from the GME (0). The probability distribution for 
the number of transferred particles is obtained from the 
n-resolved density by tracing over the system degrees of 
freedom. 



P{n,t) = TT{p{n,t)}. 



(2) 



Obviously, probability must be conserved, such that 
^jjP(n,t) — 1. In order to study the cumulants of 
P(n, t) it is convenient to introduce a cumulant gener- 
ating function (CGF) S{x, t) via the definition 



E 



(3) 



from which the cumulants ((n™)) follow as derivatives 
with respect to the counting field x X = 0, 



d{ix) 

Alternatively, one can write 

e^(x.*)=Tr{p(x,i)}, 



(4) 



(5) 
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which defines the ^-dependent density matrix 



By going to Laplace space via the transformation 



p{x,z) = / dtp{x,t) 



(7) 



Equation (Q) transforms to an algebraic equation reading 

zp{x, z) - pix, t^Q)^ W(x, z)p{x, z) + 7(x, z). (8) 

This equation can be solved formally by introducing the 
resolvent 

g(x,z) = [z-W{x.z)]-\ (9) 

and writing 

P(X, z) - Q{X, ^)[/5(x, i = 0) + 7(x, z)]. (10) 

Finally, inverting the Laplace transfoHa using the 
Broniwich integral we obtain for the CGFl3 



-1 na-\-ioo 

— / dzTr{g(x,z)[/5(x,t = 0)+7(x,^)]}e^*, 

.,(11) 

where a is larger than the real parts of all singularities of 
the integrand. 

Equation (|l^) is a powerful formal result for the CGF, 
and, as we shall see, it also leads to practical schemes for 
calculating current fluctuations. In this work, we concen- 
trate on the zero- frequency cumulants, determined by the 
long-time limit of the CGF. The case of finite-frequency 
noisepa where the inhomogeneity 7(X: z) plays an impor- 
tant role, will be considered in future works. 



B. Notational details 

Throughout this paper we will use the superoperator 
notation previously describ ed ,, ii^ g , e j^,-^^Q[, |^(j>-jused 
in a number of other works .E§HNkM3HT^J^M Us- 
ing this notation, standard linear algebra operations can 
conveniently be performed, analytically and numerically. 
Within the formalism, the memory kernel W, the resol- 
vent Q, and other operators that act linearly on density 
matrices, are referred to as superoperators and denoted 
by calligraphic characters. Conventional quantum me- 
chanical operators, like the density matrix p, acting in 
the conventional quantum mechanical Hilbert space, can 
be considered themselves to span a Hilbert space, re- 
ferred to as the superspace. The superoperators act in 
the superspace, while conventional quantum mechanical 
operators are considered as vectors using a bra(c)ket no- 
tation, i.e., V ^ \v)), where y is a conventional quantum 
mechanical operator, and \v)) is the corresponding ket in 



the superspace. Double angle brackets are used here in 
order to avoid confusion with conventional kets. In nu- 
merical calculations, bras and kets are represented by 
vectors, while superoperators are represented by matri- 
ces. The inner product between bras and kets is defined 
as {{v\u)) = Tr{VW}. Since the involved superopera- 
tors, like W and G, are not hermitian, their eigenvalues 
are generally complex. In such cases, left and right eigen- 
vectors corresponding to a particular eigenvalue are not 
related by hermitian conjugation. The left eigenvector, 
or bra, corresponding to an eigenvalue Xk is therefore 
denoted with a tilde, e.g. ((Afc|, to avoid confusion with 
the hermitian conjugate |Afe))^ of the corresponding right 
eigenvector, or ket, |Afe)). 



III. ZERO-FREQUENCY CURRENT 
CUMULANTS 

In this section we derive the recursive method for eval- 
uating the zero-frequency current cumulants. We first 
define the zero-frequency cumulants of the current as 



((/™))^ -((n™))(i) 



d d"'S{x,t) 



dt d{ixY 



(12) 

where m ~ 1,2,.... As we shall show below, the 
cumulants of the passed charge become linear in t at 
long times such that {{n™')){t) — >■ ((J™))i, and the zero- 
frequency current cumulants are thus intensive quanti- 
ties (with respect to time). Thus, in the long-time limit 
((/™ ))/((/)) = ((n'"))/((n)), and we use these two normal- 
ized quantities interchangeably throughout the paper. 

In order to find the long-time limit of the CGF, we 
consider the formal solution Eq. (pi]). The memory ker- 
nel W(x, is assumed to have a single isolated eigen- 
value Xo{x,z), which for x, = is zero, correspond- 
ing to the stationary limit of p{t), i.e., p{t) — )■ p*'*'** 
for large t. Here, p^*^^*" is the normalized solution to 
>V(x = 0, z = 0)/5*'*'^* — 0. We exclude cases, where 
the zero-eigenvalue iS|jiegenerate due to two or more un- 
coupled sub-systems. E3 In the bracket notation p'^'^* is 
denoted as |0)). The corresponding left eigenvector can 
be found be noting that the memory kernel with x = 
conserves probability for any z. This can be inferred from 
the GME in Laplace space: For normalized density matri- 
ces with Tr{p(x = 0,t)} = 1, we have Tr{p(0,z)} = 1/z, 
and Eq. (||) yields 



Tr{W(0, z)/5(0, z)} + Tr{7(0, z)} = 0. 



(13) 



It is generally possible to choose an initial state such 
that Tr{7(0,z)} = 0. The kernel does not depend on 
the choice of initial state and since Eq. dH) holds for 
any normalized density matrix p(0, z) we deduce that 
Tr{>V(0,2:) •} = 0. In the bracket notation this equality 
can be expressed as ((0|>V(0, 2;) = with the left eigen- 
vector ((0| in the superspace corresponding to the iden- 
tity operator 1 in the conventional Hilbert space. This 



5 



moreover implies that0 

Ao(0, z) — for all z. 



(14) 



We next examine the eigenvalue Ao(X: which we as- 
sume evolves adiabatically from Ao(0,0) — with small 
X and z. It is convenient to introduce the mutually or- 
thogonal projectors 



V(x, z) - V\x, = |0(x, ^)))((0(x, z)\ (15) 



and 



Q{x.z) = Q?{x.z) = l-V{x.z) 



(16) 



with V{x,z) developing adiabatically from 7^(0,0) = 
|0))((0| for small x and z. Here, ((0(x, z) \ and |0(x, z))) are 
the left and right eigenvectors corresponding to Ao(Xj z), 
which develop adiabatically from ((0| and |0)), respec- 
tively. In terms of 7'(x, z) and Q(Xi z) the memory kernel 
can be partitioned as 

W(x, z) = Ao(x, z)r{x, z) + Q(x, z)W{x, z)Q{x, z). 

(17) 

In deriving this expression we used 

P(X, z)W(x, z)V{x. z) = Ao(x, z)V{x. z). (18) 
Using the partitioning, Eq. (p^, the resolvent becomes 
V{x,z) ^, . 1 



Q{x,z) 



z- Ao(x,^) 



Q{x,z) 



z-W{x,z) 



Q{x,z) 



(19) 

For X = the first term of the resolvent has a simple 
pole at 2 = 0, which determines the long-time limit, i.e. 
it corresponds to the stationary state jf^^'^ . We denote 
the pole at z = by zo- All singularities of the second 
term have negative real parts and do not contribute in the 
long-time limit. Again, we assume adiabatic evolution of 
the pole zo(x) from zo(0) = with spail x, such that 
Zq{x) is the particular pole that solvesEZlESl 



zo - Ao(x, 2o) = 0. 



(20) 



With small x, the other singularities still have more neg- 
ative real parts and the pole zo(x) again determines the 
long-time behavior. From Eq. (O) we then find for large 
t 

eS(x,t) ^ 7^(^^^p)g2o(x)t^ (21) 

where 

i^(x,zo) -Tr{7'(x,^o)[/5(x,t = 0)+7(x,^o)]}. (22) 

From the definition of the zero-frequency current cumu- 
lants in Eq. (f2|) we then establish that 



zq 



(x) = E 



(23) 



We note that the CGF in the long-time limit and thus 
the zero- frequency cumulants do not depend on the initial 
state /5(x,i = 0) and the inhomogeneity 7(x, z). In con- 
trast, both p(x, t — Q) and 7(x, z) must be appropriately 
incorjiQrated in order to calculate the finite-frequency 
noise. Ea Equations (^ ) and (^3|) form the main theo- 
retical result of this septkua, generalizing earlier results 
for Markovian systems.cElEj In the Markovian limit, the 
memory kernel and the corresponding eigenvalue close to 
have no z-depaodence, and Eq. (|o|) immediately yields 
zoix) — ^o{x)f^^ where Ao(x) is the eigenvalue of the 
z-independent kernel, which goes to zero with x going to 
zero, i.e. Ao(0) — 0. 

Although, we have formally derived an expression for 
the CGF, it may in practice, given a specific memory 
kernel W(x, z), be difficult to determine the eigenvalue 
Ao(x, z) including its dependence on x and z. Moreover, 
the solution of Eq. (^^ itself poses an additional problem, 
which needs to be addressed in the non-Markovian case. 
In the Markovian limit, only derivatives of the eigenvalue 
Ao(x) with respect to the counting field x need to be 
determined. However, with the superoperator W(x) be- 
ing represented by a matrix of size N x N, there is no 
closed- form expression for the eigenvalue Ao(x) already 
with iV > 4. The immediate alternative strategy would 
then be to calculate numerically the eigenvalue and the 
derivatives with respect to x- Typically, however, this is 
a numerically unstable procedure, which is limited to the 
first few derivatives EJ Consequently, we devote the rest 
of this section to the development of a numerically sta- 
ble, recursive scheme that solves Eqs. ( pO| ) and (23) for 
high orders of cumulants, including in the non-Markovian 
case. 



A. Recursive scheme 

1. The Markovtan case 

We consider first the Markovian case,E3^ before 
proceeding with the general non-Markovian case. In 
the Markovian case, the memory kernel W has no z- 
dependence, and the current cumulants are determined 
by the eigenvalue Aq (x) which solves the eigenvalue prob- 
lem 



W(x)|0(x))) =Ao(x)|0(x))), 



(24) 



where Ao(0) = 0. Wc find the eigenvalue using perturba- 
tion theory in the counting field x in a spirit similar to 
that of standard Rayleigh-Schrodinger perturbation the- 
ory. To this end we introduce the unperturbed operator 
W = yV(0) and the perturbation AyV(x) such that 



W(X) = W + AW(X). 
We can then write 

Ao(x) = ((6|AW(x)|0(x))), 



(25) 



(26) 
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where we have used ((0|>V = and chosen the conven- 
tional normahzation ((0|0(x))) = 1. We moreover employ 
the shorthand notation V ^ V"^ = 7^(0, 0) = |0))((0| and 
Q = = 1 — for the projectors introduced in the 
previous section, and write 

|0(x))) = |0)) + Q|0(x))), (27) 

consistently with the choice of normalization. Using that 
W = QWQ, Eq. (H) can be written 

QWQ|0(x)» = [Ao(x) - AW(x)]|0(x))). (28) 



Next, we introduce the pseudo-inversa2jj2J defined asti2 

n^QW-^Q. (29) 

The pseudo- inverse is a well-defined object, since the in- 
version is performed in the subspace corresponding to Q, 
where W is regular. By applying TZ on both sides of Eq. 
(H) we find 



As illustrative examples we evaluate the first three cur- 
rent cumulants using the recursive scheme, 

((/^))m=((0|W(i)|0)), 

((/2)>M =((o| _ 2w(l)7^w(l)) |0)), 

((/')) A/ =((o| (w(') + 6w(i)7ew(i)7ew(i) (34) 

-6((/l))A^W(l)7^2w(l)) |0)), 

having used |0)) = |0(o))) and n\<S)) = 0, since Q|0)) = 0. 
The subscript M reminds us that these results hold for 
the Markovian case. The expressions ( p^ ) for the first 
three cumulants are equivalent to the ones derived in Ref. 

, albeit using a slightly different notation. Importantly, 
the recursive scheme presented here allows for an easy 
generation of higher order cumulants, either analytically 
or numerically. 



Q|0(x))) = 7^[Ao(x) - AW(x)]|0(x))), (30) 

which combined with Eq. (^^ yields 

|0(x))) = |0)) + 7^[Ao(x) - AW(x)]|0(x))). (31) 

Equations (|2^) and ( ^l[ ) form the basis of the recursive 
scheme developed below. 

We first Taylor expand the eigenvalue Ao(x): the eigen- 
vector |0(x))), and the perturbation AyV(x), around 
X = as 



^o(x)=f:^((n), 



oo 



|o(x)))=E 



oo 



AW(X)-E^W'"^ 



|o("))), 



(32) 



where we have used that Ao(0) — and AWfO) = 0. 
Inserting these expansions into Eqs. (^6|) and (|3^), and 
collecting terms to same order in x, we arrive at a recur- 
sive scheme reading 



((/"))= ((o|w(™)|o("-"))). 



|o(")))=7^E 



(33) 



((/")) - w^")] lo^"-"))), 



for n = 1, 2, . . .. The recursive scheme allows for system- 
atic calculations of cumulants of high orders. 



2. The non-Markovian case 

We now proceed with the non-Markovian case, where 
we first need to consider the eigenvalue problem 



W(x,^)|0(x,z))) = Ao(x,^)|0(x,^))) 



(35) 



where Ao(x, 2) is the particular eigenvalue for which 
Ao(0, z) = 0. The basic equations, Eqs. (|6|) and (31), are 
still valid, provided that Ao(x), |0(x))), W, and AW(x) 
are replaced by Ao(x,^;), |0(x, 2))), W = >V(0,0), and 
AW(x, z) = >V(x, z) - >V(0, 0), respectively, i.e.. 



Ao(x,^) = ((0|AW(x,z)|0(x,z))), 



(36) 



and 
|0(X,^))) 



-7^[Ao(x,^)-AW(x,z)]|0(x,z))). (37) 



Again, we Taylor expand all objects around x — 0, but 
in this case also around z = 0, 



AW(X,^) 



(*X) 



n 



n! l\ 



n,l=0 



n,l=0 



(38) 



E 

n,l=0 



i\ n 



with yv(o,o) ^ by definition and c^°'^'> = 0, since 
Ao(0,z) = 0. Inserting these expansions into Eqs. ( |3^ ) 
and ( |37| ) and collecting terms to same orders in x and z, 
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we find a recursive scheme reading 

I " \ I ' ^ ^^Q|yy(m,fe)|Q(n-,n,i-fc) 



m— 1 



E 

fc=0 



\m. 



I 

Eu 



In case the memory kernel has no z-dependence, corre- 
sponding to the Markovian case, only terms with I — 
are non-zero, and the recursive scheme reduces to the one 
given in Eq. (33). In particular, the coefficients c'-"'°^ 
equal the current cumulants ((/")) a/ in the Markovian 
limit of the kernel, z 0. 

In the non-Markovian case, we need to proceed with 
the solution of Eq. (|2^) for zq and extract the current 
cumulants ((/")). Inserting the expression for zq in Eq. 
( p3| ) into Eq. (po| ) and using the expansion of Ao(x, 
given in Eq. (p8[), we find 



71=1 



an 



V Mi 

^ k\ n 

k,l=0 




Collecting terms to same order in x, we find 



k\ l\ 

k,l=0 

in terms of the auxiliary quantity 

((/"0> ((/"0> 



pik,l) 



k 

E 

rii ,...,n, = l 

711 + .. .+ni=k 



I > 1, 



(40) 
(41) 

(42) 



where only terms in the sums for which ni -\- . . . -\- ni = k 
should be included. For / = 0, we have p('^'O) = 6^^- 
The auxiliary quantity can also be evaluated recursively 
by noting that 



k 

p{k,i) ^ , 

n=l 



p{k — n.l — l) 



(43) 



with the boundary conditions pC^^") = 5k,o, P(°'') = Jq,;, 
and pC^'-i) = 0. 

When combined, Eqs. (|3|, |3|) constitute a recur- 
sive scheme which allows for numerical or analytic cal- 
culations of cumulants of high orders in the general non- 
Markovian casCf— As simple examples, we show the first 
three cumulantsc3 obtained from Eqs. (^Tj), (p3|), in terms 
of the coefficients c*^"''-* 

((/3)) ^c^^fi) + ic^m^ii,!) (44) 
+ 3c(i,o) L(i.o)c(i'2)+2(c(i'i))2-fc(2.i)" . 



In general, the n'th current cumulant ((/")) contains the 
coefficients 



5(';^)9iAo(x,^)lx,.- 



(45) 



^(m,fc) _yy(m,fc) |g(n-m,i-fe)^^^ 

(39) 



with 1 < k+l < n. However, coefficients of the form c^^''-* 
are zero since Ao(0,z) = as discussed below Eq. (13) 
and it thus suffices to consider I < n — 1. From Eq. ( 39) 
it follows that c'-'^''' depend only on >V(™'") with m < k 
and n < ? so that we can conclude that the n'th cumu- 
lant of the current depends at maximum on the (n— l)'th 
time-moment of the memory kernel dtt"■~^W{x,t)■ 
^n particular, this implies that the mean current is a 
purely Markovian quantity depending only on the time- 
integrated memory kernel while the second and higher 
order cunmlants deviate from the results in the Marko- 
vian case.cJ 

The coefficients c*^"''-' can be found from Eq. (^9|). Co- 



efficients of the form c*-"''^-* only contain zeroth order 
terms in z and are, as already mentioned, equal to the 
current cumulants ((/")) a/ in the Markovian limit, i.e.. 



c("'") = ((/")) M, n=l,2,3. 



(46) 



For the other coefficients entering the expressions in Eq. 
(^J) for the first three non-Markovian current cumulants, 
we find 



.(1,1) ^/ 



.(1.2) _/ 



+ 2W(l'0)7^W(l^o)7^W(o•l) - 2W(l'l)7^W(l'°) 

_2w(l'")7^w(l■l) - w(2'0)7ew(°'i)) |o)). 

(47) 

Again, as in the Markovian case, higher order cumulants 
including the coefficients c^*^''-* are readily generated, an- 
alytically or numerically. The results presented here can 
be generalized to the statistics of several different counted 
quantities as in Ref. |9^,^, and cross-correlations can be 
evaluated pusing the same compact notation developed in 
this workEa 

B. Notes on evaluation 

As previously mentioned, the size of the memory ker- 
nel yV(x, z) could in practice hinder the calculation of 
Ao(Xj z) and the solution of Eq. (pO[), and thus the eval- 
uation of the current cumulants. The recursive scheme 
described above, however, only relies on the ability to 
solve matrix equations and perform matrix multiplica- 
tions. Both of these operations are numerically feasible 
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and stable, even when the involved matrices are of large 
dimensions. In general, the recursive scheme requires the 
following steps: The stationary state must be found by 
solving 



W|0)) = 0, 



(48) 



with the normalization requirement 
Trji^p''*^'} = 1. Secondly, the x and z derivatives 
of the memory kernel must be found 



(49) 



for {n,l) 7^ (0,0). Typically, the dependence on the 
counting field x enters matrix elements in an exponential 
function (see e. g. Refs. H, ^ ^ and examples in Sees. 
and|v|, e. g. as a factor of e*-^ , for which the deriva- 
tives with respect to x S'^e easily found analytically. The 
z-dependence of the matrix element [W{x, z)]kj can be 
written 



such that 



kj 



dt{~ty 



dt[Wix,t)]kje--\ (50) 



5(\)W(x,i)lx-o (51) 



The integration over time can be performed in a numer- 
ically stable manner for arbitrary n,0 thereby avoiding 
taking numerical derivatives with respect to z. 

Finally, matrix multiplications have to be performed. 
Here, special attention has to be paid to terms involv- 
ing the pseudo-inverse TZ, i.e. TZ\x)), where \x)) for ex- 
ample has the form W'^^O)) in the expression for the 



coefficient c^^'^'> in Eq. (47). In order to evaluate such 
expressions we introduce \y)) as the solution(s) to 



such that 



W|y)) = Q\x)), 



Q\y))=n 



(52) 



(53) 



which can be verified by applying TZ on both sides of 
Eq. dH) and using that UW ^^ QW'^QW = Q and 
TZQ — TZ. The projector Q in Eq. ( p2| ) ensures that the 
right hand side lies in the range of W, and since W is 
singular, the equation has infinitely many solutions. The 
solutions can be written 



\yo)) + c|i 



c e C, 



(54) 



where |j/o)) is a particular solution to Eq. (p2[), which can 
be found numerically. We then obtain TZ\x)) by applying 
Q to \y)) according to Eq. (p3|) and find 



TZ\x)) ^Qi\yo))+c 
since Q|0)) = 0. 



= Q\yo)), 



(55) 



In App. ^we describe a simple numerical algorithm for 
solving Eqs. (^8|) and (^). For very large dimensions of 
the involved matrices, it may be necessary to invoke mO|Ea 
advanced numerical methods to solve these equations.EEl 
Numerically, the recursive scheme is stable for very high 
orders of cumulants (up to order ^ 100), which we have 
tested on simple models. The results presented in this 
work have all been obtained using standard numerical 
methods as the one described in App. ^ 



IV. ASYMPTOTICS OF HIGH-ORDER 
CUMULANTS 



Before illustrating our methods in terms of specific 
examples, we discuss the asymptotic behavior of high- 
order cumulants. As some of us have recently shown 
certain ubiquitous, features are expected for the high- 
order cumulants. O In particular, the absolute values of 
the high-order cumulants are expected to grow factori- 
ally with the cumulant order. Moreover, the high-order 
cumulants are predicted to oscillate as functions of ba- 
sically any parameter, as well as of the cumulant order. 
This behavior was confirmed experimentally by measure- 
ments of the high-order transient cumulants of electron 
transport through a quantum dot.Ej In the experiment, 
the transient cumulants indeed grew factorially with the 
cumulant order and oscillated as functions of time (before 
reaching the long-time limit), in agreement with the gen- 
eral prediction. For completeness, we repeat here the es- 
sentials of the theory underlying these asymptotic prop- 
erties of high-order cumulants. 

The asymptotic behavior of high-order cumulants fol- 
lows from straightforward considerations. In the follow- 
ing we denote the CGF by 5(x, {A}), where {A} rep- 
resents the set of all parameters needed to specify the 
system; whether the dynamics is Markovian or non- 
Markovian is irrelevant. In general, we can assume that 
the CGF has a number of singularities in the complex-i^ 
plane at ix = iXjy 3 = 1,2,3..., which can be either 
poles or branch-points. Typically, the positions of the 
singularities depend on {A}. Exceptions, where the CGF 
has no singularities, do exist, e. g. the Poisson process, 
whose CGF is given by an exponential function, but we 
exclude such cases in the following. 

Close to a singularity ix — iXj-, we can write the CGF 

as 



^(X,{A}) 



(56) 



for some Aj and , determined by the nature of the sin- 
gularity. For example, for a finite-order pole fij denotes 
the order of the pole, while fij = — 1/2 would correspond 
to the branch point of a square-root function. Logarith- 
mic singularities can be treated on a similar footing with 
only slight modifications.Eil The derivatives with respect 
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to the counting field are now 



diixY 



with 



B. 



+ 1) ■ ■ ■ Uij+m- 1) 



(57) 



(58) 



for m > 1. As the order m is increased this approx- 
imation becomes better away from the sin. 



jdtjy at 

X — Xj according to the Darboux theorem.tSE!i!l'E23 For 
sufficiently high m, the cumulants of the passed charge 
can thus be written 



a™5(x,A) 



A B 



-i(m+A'j) arg(ixj) 



(59) 



\m+iij 



where the sum runs over all singularities of the CGF. 
Here, we have written the singularities as 



(60) 



where \ixj\ is the modulus of the singularity ixj and 
arg(zxj) is the corresponding complex argument. In gen- 
eral, the singularities ixj together with the factors Aj 
come in complex conjugate pairs, ensuring that the ex- 
pression in Eg. (p9|) is real. 

From Eq. (|59| ) we deduce that the cumulants grow fac- 
torially in magnitude with the order ni due to the factors 
Bm,iij given in Eq. (^). We also see that the high-order 
cumulants are determined primarily by the singularities 
closest to zero. Contributions from other singularities are 
suppressed with the relative distance from zero and the 
order m, and can thus be neglected for large m. Impor- 
tantly, we observe that the high-order cumulants become 
oscillatory functions of any parameter among {A} that 
changes argfix^) as well as of the cumulant order m [see 
also Eq. (61) below]. We refer to these ubiquitous fea- 
tures, which should occur in a large class of transport 
processes, as universal oscillations. For example, we ex- 
pect oscillations of high-order cumulants for basically any 
transport process described by a GME, since the CGF for 
these systems typically have logarithmic singularities at 
finite timepti or square-root branch points in the long- 
time limitJl^ Factorial growth and oscillations as func- 
tions of various parameters can be founcL-4it^Y®*|al inde- 
pendent studies of high-order cumulants,ca'E£2rt2j as well 
as in the recent experiment described in Ref. demon- 
strating the generality of the phenomenon. Similar ob- 
servations and discussions can alsp-beJound in quantum 
opticsEij and high-energy physics, Eilrti3 further confirm- 



ing the prediction. We note that in the long-time limit, 
the positions of thjs dominating singularities are no longer 
time dependentjlj and the cumulants cease to oscillate as 
functions of time. Instead, the cumulants of the passed 
charge become linear in time, as previously discussed in 



A simple (and common) situation arises if only 
two complex conjugate singularities, |ixo|e*'^''^*-^° and 
l*Xo|e~''^'^^*-^°, are closest to zero. In that case. Equa- 



tion (|59|) immediately yields 



2\Ao\B, 



cos[(m-f /xo) argixo-argAo] . (61) 



Sec. [II 



Using this expression we can determine the positions of 
the dominating singularities from numerical calculations 
of the high-order cumulants as we shall demonstrate in 
the second example considered in Sec. ^ We note that 
while the factorial growth and the oscillations are system 
independent, other features, for example the frequency of 
the oscillations, are determined by the particular details 
of the system under consideration. 

Finally, we mention the Eetmn-Frobenius theorem re- 
garding stochastic matricesL^ltiJ which implies that the 
CGFs considered in Sec. ^ must be analytical functions 
at least in a strip along the real axis in the complex-ix 
plane. This has important consequences especially for 
the nature of the high-order cumulants which rests heav- 
ily on the analytical properties of the CGF. We illustrate 
this statement in both examples in Sec. 0. 



V. MARKOVIAN SYSTEMS 

A. Electron bunching in a two-level quantum dot 

In our first example we study electron bunching in 
transport through a two-level quantum dot as described 
by Belzig in Ref. Q Due to the relatively simple ana- 
lytical structure of the model, it is possible to illustrate 
the concepts of universal oscillations introduced above. 
The model allows us to test the accuracy of our numeri- 
cal calculations of high-order cumulants against analytic 
expressions. 

We start by summarizing the setup in Belzig's model. 
Consider a single quantum dot with two single-particle 
levels coupled to voltage-biased source and drain elec- 
trodes. The two levels serve as parallel transport chan- 
nels. Due to strong Coulomb interactions on the quan- 
tum dot only one of the levels can be occupied at a time. 
The system exhibits super-Poissonian bunching trans- 
port in cases where both levels are coupled by the same 
rate F^ to the, say, left lead, whose Fermi level is kept 
well above both levels, while the couplings to the other 
lead are markedly different, such that one level is coupled 
to the right lead by the rate <^Tl and the other by 
xV n with x <C 1. This situation can arise for example, 
if the two levels are situated above and slightly below, 
respectively, the Fermi level of the right lead at a finite 
electron temperature. 

This particular configuration leads to bunching of elec- 
trons in the transport due to the existence of the blocking 
state: if the dot is empty there is equal probability for 
either of the two levels to be filled. Current runs easily 
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through the first level, while the other level effectively 
is blocked, or niore precisely, the transport through the 
level is limited by the very small right rate xTji, con- 
stituting a bottleneck. The transport thus proceeds in 
bunches of electrons passing intermittently through the 
first level separated by quiet periods of blocked transport 
when the other level is occupied. This bunching effects 
leads to super-Poissonian noise with a Fano factor above 
unity. For more detailed discussions of the model as well 
as its generalizations to many levels, the reader is referred 
to Ref. ||. 

The counting statistics of the system can be obtained 
from a Markovian rate equation for the probability vec- 
tor p = (poiP+iP-)"'", containing the (n- resolved) prob- 
abilities po,-,+ for the quantum dot to be empty, or the 
first non-blocking) or second (— , blocking) level being 
occupied, respectively. The corresponding x-dependent 
rate matrix reads 



W(x) 




(62) 



Here, we have rescaled the time and set F^ = 1 while re- 
naming Tfi = F in order to simplify the analytic results 
in the following. We have also made a minor modification 
of the model in Ref. ^ by including the back-flow into 
the blocking level from the right lead. This modification, 
however, changes only slightly the detailed quantitative 
results, while leaving the main qualitative features iden- 
tical in the limit of interest a;, F <C 1. 

Since the model involves only three states, the CGF 
can be found analytically in the long-time limit. The 
full expression is too lengthy to be presented here, but 
in the limit x, F ^ 1, it reduces to the result by BelzigE2 
(also for our slightly modified model; note, however, the 
opposite sign convention for the CGF in Ref. BSh 



1 



2 - e»x 

Clearly, the CGF has simple poles at 

iXj ^\n2+ j2TTi, j = ..., -1,0,1, . 



(63) 



(64) 



with the pole ixo = ln2 being closest to 0. However, 
according to the Perron- Frobenius theorem mentioned in 
Sec. ^ the CGF cannot have singularities on the real 
I X-axis. 

In order to illustrate this point, we consider the ex- 
pected behavior of the high-order cumulants based on 
the CGF above. Close to the singularity ixo, we approx- 
imate the CGF by the first non-zero term of the Laurent 
series 



six,t) 



Txt 



(65) 



1X0 - «X 

This corresponds to Eq. ( |5^ ) with Aq — Txt and fio — 
1. From Eq. (^9|) we then obtain a simple asymptotic 
expression for the high-order cumulants reading 



{{n)/rx 


m = 1 


2 


3 


4 


5 


6 


Single-pole approx. 


2.000 


6.000 


26.00 


150.0 


1082 


9366 


Single-pole asympt. 


2.081 


6.006 


25.99 


150.0 


1082 


9366 


Numerics 


1.978 


5.880 


25.18 


143.3 


1017 


8644 



TABLE I: Normalized zero- frequency current cumulants for 
transport through a two-level quantum dot. Single-pole ap- 
proximation results have been obtained by direct differentia- 
tion of the CGF in Eq. (^3]) or its asymptotic expression Eq. 
(^6|), respectively. The numerically exact results have been 
obtained using our recursive scheme and the rate matrbc in 
Eq. (S) with X = 0.001 and T = 0.01. 



Here, the subscript indicates that the expression has 
been obtained using the approximate CGF in Eq. ( |6^ ) 
with only a single singularity closest to zero. In Table 
Q we compare the asymptotic expression with results for 
the first six cumulants obtained by direct differentiation 
of the CGF in Eq. (|63|). The asymptotic results are very 
close to the exact derivatives of the approximate CGF. 
Despite the good agreement with the approximate re- 
sults, the asymptotic expression in Eq. (|66| ) does not re- 
produce our numerically exact results, also shown in the 
table, obtained using our recursive scheme. In particu- 
lar for high orders, the asymptotic expression starts to 
deviate significantly from the numerically exact results. 

As anticipated above, these deviations can be traced 
back to the expression in Eq. (|6^), that we obtained in 
the limit x,F ^ 1. In order to proceed from here, we 
return to the full expression for the CGF in the long-time 
limit (not shown). A careful analysis reveals that, in fact, 
there is a pair of complex conjugate singularities closest 
to zero, and not just a single pole. The two singularities, 
denoted as ixo and (iXo)*, correspond to branch points 
of a square-root, and for small a; <C 1 the position of the 
branch point ixo is 



^Xo = ln(2 + r)-2.i±iM+4.V^.l + ^ 



(2- 



r 



(67) 



{{in)is/rx = ((n"))i,/Fxt ~ m!/(ln2 



,m+l 



(66) 



Clearly, for small a;, F <C 1 the branch points are close 
to the position of the single pole ixo = ln2. However, 
for any finite x, the two branch points have small, but 
finite, imaginary parts thus complying with the Perron- 
Frobenius theorem. The singularity structure around 
the branch point ixo is characterized by Eq. ( ^6| ) with 
/io — —1/2 and Aq « Ft-^e*'^/'*, and we can then use 
the asymptotic expressions in Eq. ( |6l[ ) for the high-order 
cumulants. In the left and central panels of Fig. || we 
compare this expression, and the single-pole approxima- 
tion in Eq. (^6|) , with numerically exact results obtained 
using our recursive scheme for F — 0.01 and two different 
values of x = 0.001, 0.01. 

Figure ^ shows several important features. Firstly, the 
(scaled) high-order cumulants indeed behave in an os- 
cillatory manner as function of the cumulant order m, 
which coincides with the cosine part of Eq. ( |6l] ) . Obvi- 
ously, for smaller x = 0.001 the period of the oscillations. 
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FIG. 2: High-order cumulants and large deviation function for bunching transport through a two-level quantum dot. Left and 
central panels show comparisons between exact numerics and the single pole approximation stemming from Eq. (^) for two 
different values of a; = 0.001, 0.01 and F = 0.01. The asymptotic expression in Eq. ( |6l| ) based on a pair of complex conjugate 



singularities is shown with full lines. Notice that B„ 



-1/2 



< 0. The right panel shows a comparison of the large deviation 



function (LDF) obtained from exact numerics and the single pole approximation in Eq. (vW, respectively. 



determined by argixo, is longer in accordance with Eq. 
(|6l|). Furthermore, for the small value of a: = 0.001, the 
asymptotic form of the high-order cumulants is reached 
around m ~ 30, while the single-pole approximation 
agrees well for lower orders, to < 10. For the higher value 
of a; = 0.01, significant deviations from the single-pole be- 
havior begin already for the fourth cuniulant, while the 
asymptotic oscillatory form holds from around to = 12. 
Notice the importance of the exact analytical knowledge 
of the singularities — even though x = 0.01 ^ 1 (to- 
gether with F = 0.01 <^ 1) may seem a very small num- 
ber justifying the usage of the single pole approximation, 
we see from Eq. that the imaginary part of the pole 
and its argument scale like ^/x = 0.1, thus invalidating 
the single-pole approximation far earlier than expected 
from a linear-in-x scaling assumption. 

A complementary view on the charge transport static 
tics is provided by the large deviation function (LDF),Ej 
which quantifies deviations of measurable currents from 
the average value. The LDF is obtained from the proba- 
bility distribution 



1 

2^ 



dxe 



S(x,t)-inx 



(68) 



and is defined as the long-time limit of \n[P{I,t)]/t, 
where / = n/t is the current. For long times, we have 
S{x, t) — )■ Ao(x)^ and the integral can be evaluated in the 
saddle-point approximation with the saddle-point x — Xo 
given by the solution to the saddle-point equation 



(69) 



The saddle-point equation implies a parametric depen- 
dence of the saddle-point xo = Xo(-^) on the current /. 
Using the saddle-point approximation, the LDF becomes 



ln[P(/,t)] 
t 



Ao(xo) - ilXo- 



(70) 



We first solve the saddle-point equation for the approxi- 
mate CGF in Eq. iM) and find 



ln[Pi,(/,t)] VI + 8k -3 



{I)t 



Klog 



16k 



(1 + ^/TT8^)2 

(71) 

where k — //(/) and the subscript again reminds us 
that the expression has been obtained using the approx- 
imate CGF with only a single singularity closest to zero. 
Obviously, the current must be positive {k > 0), since 
transport is unidirectional. 

Also for the LDF, we can compare the analytic approx- 
imation with numerical exact results. To this end, we 
need to solve the saddle-point equation taking as start- 
ing point the kernel in Eq. (^2|). The derivative of the 
eigenvalue Ao(x) is now calculated using the Hellman- 
Feynman theorem, writing 



KU) = |^((o(x)|w(x)|o(x))) 
= ((o(x)|w'(x)|o(x))>, 



(72) 



where ((0(x)| and |0(x))) are left and right eigenvectors 
of W{x}, respectively, corresponding to the eigenvalue 
Ao(x), and ((6(x)|0(x))) = 1- For a given value of x 
we calculate numerically the left and right eigenvectors 
((O(x)l and |0(x))) and find Aq(x) using the expression for 
the derivative in Eq. ([72|). With this procedure we search 
numerically for the value of x = Xo that solves Eq. ( |69| ) 
for a given value of /, and with the solution xo we eval- 
uate the LDF using Eq. (|70|). We find that xo is purely 



12 



imaginary.c2l We note that, in principle, the existence of 
a saddle point solution is not guaranteei-jii, the whole 
range of currents, and there are examples |22I'E3 where the 
behavior of the LDF changes abruptly at finite values of 
/ due to singularities of the CGF on the real i^-axis. In 
our case, however, the Pcrron-Frobenius theorem ensures 
that the CGF is analytic on the real z^-axis and the LDF 
is smooth as function of /. In the right panel of Fig. ^ 
we show a comparison between exact numerics and the 
analytic result (|7l| ) for the large deviation function in 
the single pole approximation. Around the mean value 
/ ~ (/) the analytic result agrees well with numerics. 
However, in the tails of the distribution a clear disagree- 
ment between the analytic approximation and numerics 
is visible. The disagreement reflects the deviations for the 
cumulants seen in the central panel of Fig. ^ We remark 
that measurements of the LDF recently have become af^ 
cessible in experiments on real-time electron counting.tS 
The discussion in this subsection illustrates the need 
for careful considerations when manipulating CGFs an- 
alytically. Concerning cumulants, we deal with two op- 
posite and non-commutative orders of limits: for a fixed 
order of cumulants, a limiting procedure with changing 
parameters converges to the approximate form given by 
the appropriate limit of the CGF, such as the single-pole 
approximation in Eq. (|63| ) in our case. However, the 
convergence of the CGF is not uniform in x due to po- 
tential singularities and thus for fixed parameters, high 
order cumulants generically take on the universal oscil- 
latory form discussed above. One should thus be careful 
when using limiting forms of a CGF to extract cumulants 
of arbitrary orders. In general, the low-order cumulants 
follow the predicted pattern reasonably well, but at some 
point significant deviations appear and the universal os- 
cillatory behavior should emerge. The order at which 
this crossover occurs depends on details of the analytical 
structure of the CGF and may be hard to predict. As we 
have shown explicitly, deviations of the cumulants from 
exact results are also clearly visible in the large deviation 
function. 



B. Transport through a vibrating molecule 

In our next example, we consider a model of charge 
transport _tiux)nflh_|a, molecule coupled to quantized 
vibrations £30'c3c3l3 In the regime of weak coupling to 
the electronic leads, electron tunneling can be described 
using Fermi's golden rule rates for transitions between 
different vibrational and charge occupation states. For 
strong electron-phonon coupling, the large shift of the 
oscillator equilibrium position due to an electron tunnel- 
ing onto the molecule suppresses the (Franck-Condon) 
overlap between the initial and final vibrational state for 
low-lying oscillator states. This leads to surpressed tun- 
nel rates at low bias-voltages, so-called Franck-Condon 
blockade. For larger bias- voltages, higher-excited oscilla- 
tor states become available, and the system can escape 
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FIG. 3: (color online). Transport through vibrating molecule. 
The molecule is coupled to the left (right) lead with coupling 
Fl (Tr)- The bias difference eV = I-Il — f^R drives single 
electrons through the molecule. The system is operated in 
the Coulomb blockade regime, where only m = or m = 
1 additional electrons are allowed on the molecule. As an 
electron tunnels onto the molecule, the equilibrium position 
of the molecule is shifted due to the electric field E. The two 
harmonic potentials corresponding to m = 0, 1 are shown. 
The damping rate of the vibrating molecule is denoted as K. 



the blockade regime. For weak oscillator dampings, sev- 
eral electrons can be transferred through the molecule, 
once the blockade is lifted, until a charge transfer event 
eventually leaves the oscillator in the ground state and 
the current is suppressed again. Such dynamical Franck- 
Condon blockade processes have been predicted to lead, 
to very large enhancements of the zero- frequency noise.cil 
Recently, Frank-Condon blockade was observed in expjet. 
iments on suspended carbon nanotube quantum dots.ti^ 

The system considered in the following is depicted in 
Fig. |. Here we follow to a large extent the description 
of the model given in Refs. The Hamiltonian of 

the system and the detailed derivation of the resulting 
Markovian GME are given in App. where the vari- 
ous parameters of the model are also defined. Due to 
the large number of oscillator states, there is little hope 
for obtaining a closed-form expression for the CGF that 
would allow for any analytic manipulations. Instead, as 
we shall see, the numerically evaluated high-order cu- 
mulants can be used to extract the precise location of 
the dominating singularities of the CGF. We concentrate 
in the following on the unequilibrated oscillator regime, 
where the damping rate of the oscillator is much smaller 



than the electron tunneling rates, K ^ ^l/r- As ex- 
plained above, the combination of strong electron-phonon 
coupling and weak oscillator damping leads to dynamical 
Franck-Condon blockade, resulting in a large enhance- 
ment of the current noise as demonstrated in Ref. ^ 
using Monte-Carlo simulations. In Ref. ^ the analysis 
was extended to the full distribution of the transferred 
charge and an analytic approximation for the CGF was 
presented based on an avalanche-type of transport, where 
"quiet" periods of transport are interrupted by a sequel 
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of self-similar charge avalanches. The analytic result for 
the CGF was shown to agree very well with Monte-Carlo 
simulations of the probability distribution P{n, t). How- 
ever, similarly to the previous example, the approximate 
CGF has a single, simple pole on the real-«x axis, vio- 
lating the required properties of the CGF, mentioned at 
the end of Sec. LV, thus making it unsuited for predic- 
tions of the high-order cumulants. In particular, within 
this approximation, the high-order cumulants would not 
oscillate, which contradicts our numerical findings. 

Oscillations of the high-order cumulants with system 
parameters must be due to singularities located away 
from the real-zx. In the following, we assume that 
the CGF has a pair of complex-conjugate singularities, 
«Xo = |«Xo|e''"'8*'^° and Izxole"""^*-^", closest to zero. As 
we will now show, the positions of these singularities can 
be found from our numerical calculations of high-order 
cumulants. To this end, we define 



ao = Ao/iixor" 
and rewrite Eq. (|6l| ) as 

//^ni\\ ^ 2|ao|i3„i,p(, 



cos[77^arg^xo-arg ooj 



(73) 



(74) 



Following the ideas of Rcf. 11£, Sec. 4, we find for the 
ratios of two successive cumulants 



((n-)) 



and 



- cos[argzxoj 

sin[arg ixo] tan[m arg ixo - arg ao] 

(75) 



: cos[argixoj 

sin[argixo] tan[margixo ~argao] . 

(76) 



Adding the two left and right hand sides, respectively, 
and rearranging, we obtain the equation 

2(m + /.o)((n"))Nxo| cos[arg»xo] - |^xo|' = 

((n™-i))(m + /io-l)(m + Mo). 

(77) 

Using the substitution m — )• m + 1, we obtain an ad- 
ditional equation and thus arrive at a linear system of 
two equations which we solve for |zxo| cos[arg ixo] and 
|zXoP and thereby find ixo- The method takes as in- 
put ((n™)), and and the ac- 
curacy is-escpected to improve with increasing cumulant 
order m.ti3 

Having determined ixo, we find ao in a similar spirit 
by rewriting Eq. (74) as 



.2B,„,^„ [Re{(zxo)-™}Re{ao} 
-Im{(iXo)-'"}Im{ao}] . 




(78) 



FIG. 4: (color online). High-order (normalized) cumulants 
for unequilibrated molecule. Numerically exact results are 
shown together with the asymptotics described by Eq. (|6l|). 
Parameters entering Eq. (pll) are /io = —1/2, Ao — 1.4810 x 
10-^e-'''■■'^■'^ and ixo = O.OllSe'"-^^'^^ System parameters 
(defined in App. |^ are given in units of the natural oscillator 
frequency (with e, fi, fes = 1) = 3uuo, T — Tl ~ Tr = 
0.001a;o, T = 0.05a;o, K = 10~^°a;o, e = 16a;o, a = 4, C2 = 0. 
In the numerical calculations we have used = 15 oscillator 
states. 



Again, we obtain via the substitution to —> to -I- 1 a linear 
system of two equations that we solve for Re{ao} and 
Im{ao} and thus find ao = Re{ao}-l-«Im{ao}. Finally, we 
determine Aq from Eq. ([73|). More advance methods fo 
extracting the positions of singularities are available,ti^ 
but they require solutions of non-linear equations and 
will not be considered here. 

In order to extract ixo and Aq from the high-order cu- 
mulants, we need to know the nature of the singularities 
and hence /io. Typically, t he si ngularities are square- 
root branch points (see Ref. 101, Sec. 7.5) and we thus 
take /io = — 1/2. In Fig. ^ we show numerical results 
for the (normalized) cumulants as function of the order 
TO together with the asymptotic expression Eq. (^ij) for 
the high-order cumulants with ixo and Aq found using 
the method described above. The asymptotic expression 
shows excellent agreement with the numerically exact re- 
sults. For TO > 5, we see trigonometric oscillations whose 
frequency is determined by arg?xo- We note that a good 
agreement between our numerical results and the asymp- 
totic expression could only be obtained with fiQ = —1/2, 
thus confirming that the singularities stem from square- 
root branch points. 

The large deviation function can also be evaluated nu- 
merically using the method described in the previous sub- 
section. In Fig. 1^ we show numerical results for the large 
deviation function with different values of the damping 
K. For large dampings, the oscillator is essentially equi- 
librated and the measurable currents are closely centered 
around the mean current ((/)). As the damping is low- 
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FIG. 5: (color online). Large deviation function for the vi- 
brating molecule. Results are shown for different values of 
the damping K, going from the unequilibrated regime, at low 
K <^ r, to K T, where the molecule equilibrates between 
each tunneling event. For the unequilibrated much 
larger range of currents are probable, compared to transport 
through the equilibrated molecule. System parameters (de- 
fined in App. p|) are given in units of the natural oscillator 
frequency (with e,h,kB = 1) V = Sloq, Tl ~ ~ O.OOIcjq, 
T = O.OStJo, e = 16tJo, ci = 4, C2 = 0. 



ered, we approach the unequilibrated regime, where the 
transport statistics is dominated by avalanche transport 
with a corresponding large zero-frequency noise. Accord- 
ingly, the large deviation function is considerably broad- 
ened and a much wider range of currents become mea- 
surable. 



VI. NON-MARKOVIAN SYSTEMS 
A. Dissipative double quantum dot 

In the previous two examples, we focused on the 
asymptotic behavior of the high-order cumulants for two 
Markovian systems. We now turn our attention to a 
model for which a weak coupling prescription does not 
suffice and non-Markovian effects become significant. We 
focus here on the influence of memory effects on the first 
few cumulants, while referring the reader to Ref. for 
a discussion of the high-order cumulants for the non- 
Markovian system presented in this example. 

We consider a model of charge transport through a 
double quantum dot (DQD) coupled to a heat bath which 
causes dephasing and relaxation^ Such systems were 
studied experimentally in Refs. 117,118. The counting 



statistics in the transition between coherent and sequen- 
tial tunneling through DQDs has Jaeen studied theoret- 
ically by Kiefilich and co- workers. Eij In their work, de- 
coherence was described using either a charge detector 
model or via phenomenological voltage probes.Q More 
elaborate descriptions of decoherence caused by a weakly 




FIG. 6: Dissipative double quantum dot. The Coulomb block- 
aded double quantum dot consists of the left and right levels 
jL) and |_R), coherently coupled with tunnel coupling Tc and 
dealigned by e. A large bias across the system drives electrons 
through the double quantum dot from the left lead with rate 
Fl to the right lead with rate Tr. The system is coupled 
with dissipation strength q to a heat bath at temperature 
T and with Ohmic spectral function Jn{oj). The probability 
distribution of the number of transferred charges n is denoted 
P{n,t). 



coupled heat bath were given in Refs. 35 120 and shown 
to agree well with experiments. 

Here, we take these ideas further and go beyond the 
perturbative treatment of the heat bath. This situation 
has previously been investigated by Aguado and Brandes 
using a polaron transformation, assuming weak,p3upling 
to the electronic leads in the high-bias limit .E30'tHil In 
the following, we apply an alternative non-perturbative 
scheme for the coupling to the heat bath, enabling us 
to fully include broadening due to the electronic leads. 
Within this approach, we can study the cross-over be- 
tween weak and strong couplings to the heat bath and 
evaluate the effects of strong decoherence on the charge 
transport statistics. In particular, we show that only in 
the limit of weak coupling and high temperatures, the de- 
phasing caused by the heat bath can be accounted for by 
a charge detector model with a single effective dephasing 
rate. 

The model ot-ciiarge transport through a Coulomb 
blockaded DQDEiki is illustrated in Fig. |. The DQD 
is coupled to source and drain electrodes, while dissipa- 
tion is provided by an external heat bath. The DQD 
is operated in the Coulomb blockade regime close to a 
charge degeneracy point, where only a single additional 
electron is allowed on the double dot. Again, we consider 
for simplicity spinless electrons. The Hamiltonian of the 
double dot can be written 



£0 



where the pseudo-spin operators are 
^ \L){L\ - \R){R\ 

and 

= \L){R\ + \R){Ll 



(79) 



(80) 



(81) 
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respectively. Here, the two quantum dot levels \L) and 
\R) are dealigned by e and their tunnel coupling is T^- 
The energy of the 'empty' state |0) is cq. The pseudo- 
spin interacts with an external heat bath consisting of 
harmonic oscillators, 



(82) 



whose positions are coupled to the z-component of the 
pseudo-spin, adding the term VbSz to the full Hamilto- 
nian with 



(83) 



and are assumed energy- independent, such that = 
Ta, a — L, R. We count the number of electrons that 
have been collected in the right lead, and consistently 
with this choice, the counting field x has been introduced 
in the off-diagonal element of the memory kernel that 
contains the rate r^. 

The expressions for the bath-assisted hopping rates are 
derived in App. |^ and for real z they read 



(89) 



where z± = z ± is + Tf(/2. These expression are valid 
to the lowest order in the tunnel coupling Tc. The bath- 
correlation functions in Laplace space are 



Finally, the spin-boson system is tunnel-coupled to left 
(i) and right (i?) leads via the tunnel-Hamiltonian 



(84) 



ka,a=L.R 



with both leads described as non-interacting fermions, 
i.e.. 



Ha = ^^£k^c\^Ck^, a = L,R, 



(85) 



Jo 



(90) 



with 
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W{t) 
and 



doj—^{[l~cos{ujt)] coth {Puj/2)+i sin(a;t)}. 



(91) 
(92) 



kept at chemical potentials /Iq, a — L,R, and tempera- 
ture T. The full Hamiltonian then reads 



H = Hs + Ht + Hl+Hr + Hb + VbS, 



(86) 



As previously pointed out j£j'Lj the model can be mapped 
onto that of transport through a superconducting single 
electron transistor, when the charging energy is much 
larger than the Josephson energy. Throughout this ex- 
ample we take ft = A;^ = e = 1. 

As explained in App. transport through the double 
dot can be described using a non-Markovian equation of 
motion of the form in Eq. (|l|) for the three electronic 
occupations of the double dot collected in the vector 
p ~ {po, PL, Pr)^ ■ The occupation probabilities of the 
empty, left, and right states, are denoted po, pL, and 
PR, respectively. The corresponding memory kernel in 
Laplace space reads 



W(x,^) = 



Tl 






.(+) 

B 
.(+) 



rL-'(^) 



(87) 



We note that the kernel with x = has a single zero- 
eigenvalue Ao(0, z) = for all z, in agreement with Eq. 
(p^). The kernel has been derived under the assump- 
tion that the symmetrically applied bias eV =\p,L — 
between the electronic leads is much larger than the tun- 
neling rates to the leads and the temperature T. The 
tunneling rates are defined as 



r„(e) = 27r^|ifej2,5(e-£fcj, a = L,R, 



(88) 



being the spectral function of the heat bath. In this 
work we consider Ohmic dissipation characterized by a 
coupling strength a such that the spectral function reads 



Jn{oj) = 2awe-'^/'' 



(93) 



where ujc is the frequency cut-off, assumed to be the high- 
est energy scale of the system. 

In Fig. we show results for weak couplings to the heat 
bath, a <^ 1. As the two quantum dot levels are tuned 
into resonance (e = 0), the current reaches a maximum 
with a width mainly determined by T r. The correspond- 
ing values for the second and third cumulant, normal- 
ized with respect to the current, are suppressed below 
unity. The suppression is stronger for the third cumu- 
lant. Away from resonance, the mean current falls off, 
and the second and third cumulants approach unity, cor- 
responding to a Poisson process. Without coupling to the 
heat bath, a = (dotted blue line), the only broadening 
mechanism is the escape of electrons through the right 
barrier at rate Tr. This rate also defines the relevant 
energy scale to which we compare the temperature of the 
heat bath. Away from resonance, the uncoupled case 
captures well the results obtained at low temperatures, 
T < Ffl. (dashed-dotted green line). At higher tempera- 
tures, T 3> Tr (full red line), the peak in the current and 
the dips in the second and third cumulants are consid- 
erably broadened due to the strong temperature induced 
dephasing. Further results for the weak coupling limit 
are presented in Ref. [123| . 

In order to understand the behavior at high tempera- 
tures (full red line), we imagine replacing the heat bath 
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FIG. 7: (color online). Cumulants for weakly coupled heat 
bath. The first three cumulants are shown as functions of 
the dealignment e. We show results for a low (T* = 0.8r_R) 
and a high temperature (T = Sr^) as well as for the un- 
coupled case (a = 0). For the high temperature case, we 
compare with results obtained using a charge detector de- 
phasing model (short-dashed black line) and results without 
inclusion of memory effects (long-dashed black line), see text. 
Parameters are Vl = O.IFh, Tc = O.lFfl, a = 0.01, and 
ujc = ^>< lO^Ffi, and thus Fd ~ O.SFjj according to Eq. (|^. 
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by a charge detector which measures the position 



ll£,L2^ 



L25 



irons on the DQD, thereby causing dephasing 
The effects of the charge detector can be described by 
a single dephasing rate T^, entering as an additional ex- 
ponential decay of the off-diagonal elements between the 
left and right quantum dot states. As we show in App. 

this picture follows from the high-temperature limit of 
the kernel in Eq. (^^, and the corresponding dephasing 
rate is 



Frf = 2a7rT. 



(94) 



The dynamics of the system effectively becomes Marko- 
vian at high temperatures T 3> F^, where the charac- 
teristic memory time ~ {T jil2 + Tci)~^ of the kernel is 



shorter than the timescale ~ F^ over which the pop- 
ulations of the DQD evolve. In Fig. Q we see that the 
counting statistics at high temperatures (full red line) 
are well approximated by the charge detector model 
(short-dashed black line), which captures the broaden- 
ing of the peak in the current and the dips in the second 
and third cumulants. For high temperatures (full red 
line) , the large value of the dephasing rate indicates that 
the system is strongly dephased. The charge detector 
model, however, cannot account for the weak asymme- 
try between the phonon emission (e > 0) and absorption 
(e < 0) sides at low temperatures (dashed-dotted green 
line). 

In our description of the DQD system we have traced 
out the electronic off-diagonal elements, the coherencies, 
together with the electronic leads and the heat bath. Our 
derivation allows us to combine strong coupling to the 
heat bath with broadening of the electronic levels due 
to the electrodes. However, even without coupling to the 
heat bath, the kernel must be time-dependent in order to 
account for the coherent oscillations between the left and 
right quantum dot states. These coherent effects are sup- 
pressed, when the dephasing is strong, and in that limit 
we thus expect that a Markovian description would suf- 
fice. We check this assumption by plotting in Fig. ^ only 
the Markovian parts of the cumulants at high tempera- 
tures (long-dashed black line). Away from resonance, the 
Markovian parts agree well for the second and third cu- 
mulants showing that the system at high temperatures 
effectively is Markovian. The mean current, as previ- 
ously mentioned, is already a Markovian quantity andi± 
coincides with the Markovian contribution as expected.EJ 
Closer to resonance, some deviations for the second and 
third cumulants are seen as the system is not completely 
dephased. 

In Fig. 1^ we show results for larger values of the cou- 
pling to the heat bath. Also in this case, the peak in 
the current and the dips in the second and third cu- 
mulants are suppressed as the coupling is increased and 
dephasing becomes stronger. Contrary to the weak cou- 
pling regime, however, no broadening of these features 
are observed. This is not consistent with the charge de- 
tector model, which would predict an increased width 
together with suppression of the height of the current 
peak and the depth of the dips in the second and third 
cumulants. Additionally, as the coupling a is increased, 
the emission/absorption asymmetry becomes stronger as 
exchange of energy quanta with the heat bath becomes 
increasingly important. As noted before, neither the ab- 
sence of broadening nor the asymmetry of the peaks and 
dips can be accounted for by the charge detector model. 
For large values of the coupling to the heat bath, the 
DQD system completely dephases and the Markovian 
contribution describes well the behavior of the counting 
statistics, in particular away from resonance, where co- 
herent effects are less relevant. At even larger couplings, 
the heat bath tends to localize electrons to one of the two 
quantum dots, and the effective tunnel rate between the 
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FIG. 8: (color online). Cumulants in the intermediate regime 
between weak and strong coupling to the heat bath. The first 
three cumulants are shown as functions of the dealignment 
£. The values of the coupling to the heat bath are a = 0, 
0.01, 0.05, and 0.1. For large couplings (a = 0.5, 1), we com- 
pare with results obtained without inclusion of memory ef- 
fects (dot-dashed lines), see text. Parameters are Tl = r_R, 
= O.lFfl, T = 0.5Ffl, and = 5 x 10*Ffl. 



two quantum dots becomes highly suppressed. In that 
case, the current through the system is very low and the 
statistics is Poissonian. 



B. Discussion of non-Markovian systems 
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i^arious possible 
and, in partic- 



We close this section by pointing ax 
subtleties associated jwith MarkovianEZ 
ular, non-Markoviar£Hl GMEs. As we shall argue, spe- 
cial attention should be paid to the sometimes paradoxi- 
cal nature of heuristically derived non-Markovian GMEs 
(see e. g. Ref. 127) in order to ensure physically meaning- 
ful results. We moreover discuss the interpretation of the 
"mean memory time" for non-Markovian systems and we 
show that it in certain cases can turn negative. As an il- 



lustrative example, we consider unidirectional transport 
through a single electronic level. This generic model pro- 
vides us with a unifying explanation of previous results 
obtained for several different systems and, in addition, 
it displays a few possible peculiarities of non-Markovian 
transport. 

The Markovian master equation for the simple two- 
state model of unidirectional transport through^ 
electronic level is determined by the rate matris 



w(x) 




(95) 



with the counting field x corresponding to tunneling 
across the right barrier. An intuitive way to generalize 
the rate matrix to the non-Markovian case would be sim- 
ply to replace the rates F^/j^ by time-dependent rates, 
so that the rate matrix in Laplace space instead reads 



W(x,2)- 




(96) 



The non-Markovian character could be caused by exter- 
nal degrees of freedom that have been traced out, for 
example a harmonic oscillator mode coupled to th e oc- 
cupation of the electronic level, as studied in Sec. VB . 
Alternatively, it could be due to an energy-dependent 
tunneling density of states as in Ref. or many-body- 
induced effeGts-|as-in the Fermi-edge-singularity problem 
in transport ££3nH3 

The microscopic origin of the memory, however, is not 
important in the following and its effect on the current 
noise is qualitatively captured by the "mean memory 
times" 



TL,fl = - flog Fi,fl(z) I 



(97) 



If F(i) cx e"*/"^ the mean memory time has a direct phys- 
ical interpretation as the characteristic memory time r. 
However, the following statements are valid for any mem- 
ory kernel as long as it has a finite mean memory time 
T. Using Eqs. ( ^ , ^7| ) we find the followin g exp ression for 
the Fano factorjsee also Eq. (13) of Ref. 
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F = 



F2 



- 2 



(98) 



where Tl r = rL^ji{z = 0) > is the Markovian limit of 
the rates. The second term is a non-Markovian correc- 
tion to the well-known expression for the Eano factor of 
transport through a single electronic leveLu For certain 
parameter values, however, this non-Markovian correc- 
tion can make the Fano factor turn negative — a clearly 
unphysical result (the zero-frequency noise must be posi- 
tive, see e.g. Refs. ^,|o|). We discuss this issue in further 
detail below. 

Obviously, different baths producing the same mean 
memory times give rise to the same electronic noise. In- 
tuitively, one would also expect the mean memory times 
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to be positive. Since the Markovian limit of the F's 
must be positive (being rates) the non-Markovian cor- 
rection appears negative; thus the general effect of mem- 
ory on transport through a single level is a decrease 
of the noise compared to the corresponding Markovian 
limit. This statement is in line with a number of previ- 
ous findings: It explains the anomalous suppression of the 
Fano factor below 1/2 in transport through a single elec- 
tronic level coupled to a mechanical resonator reported 
in Ref. M. Additionally, the non-Markovian correction 



due to strong 
larity probleaa 
discrepancyt^ 
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atures in the Fermi edge singu- 
is responsible for the observed 
between the measured Fano factor and 
the expected result based on the Markovian part only. A 
more detailed account of this problem will be presented 
elsewhere. Finally, the suppression of the Fano factor 
compared to the Markovian approximation is also con- 
firmed by the study of an exactly solvable caseE^ in the 
regime where the non-Markovian GME provides a good 
approximation to the exact dynamics. 

Although the above expression ( |98| ) provides a uni- 
fying explanation of these three examples, it obviously 
cannot be correct in general as mentioned above. For 
weak non-Markovian behavior, where the t's are small, 
the Fano factor stays positive, and the non-Markovian 
corrections lead to a reduction of noise. In general, how- 
ever, there is no guaranty that the Fano factor in Eq. 
( p8| ) is always non-negative. This can be traced back 
to the heuristic inclusion of the non-Markovian kernel 
in Eq. (p^). While the non-Markovian kernel for uni- 
directional transport through a single level in general 
may be written in the form ( p6[ ) without counting fields, 
the inclusion of the counting field must be carried out 
carefully, for example by using well-controlled system- 
atic derivation jmpcedures, such as those based on pertur- 
bation theoriesH] Non-Markovian GMEs may, however, 
still lead to unphysical results, when employed outside 
their regime of rvalidity, as recently discussed by Zedler 
and co-workersH In the Markovian case, the heuristic 
addition of the counting field in W usually leads to cor- 
rect results fop-t|hc counting statistics, although excep- 
tions do exist In the n on-Markovian double dot sys- 
tem studied in Sec. VIA, the counting field enters the 



Markovian (z-independent) part of kernel in Eq. (87), 
and the inclusion of the counting field does not lead to 
any of the issues discussed above. 

Finally, we discuss another subtlety associated with 
strongly non-Markovian systems. Under certain circum- 
stances the mean memory time r, defined in Eq. ( p7| ) 
may in fact become negative. This happens for exam- 
ple f or th e dissipative double quantum dot studied in 
Sec. VIA. For a sufficiently small dissipation rate, we 
find for the bath-assisted rates dzT^^\z = 0) > 0, result- 
ing in negative mean memory times. Formally, there is 
no problem associated with this phenomenon (the GME 
still describes a positivity-preserving evolution), but the 
physical interpretation of the non-Markovian corrections 
is less clear due to this counterintuitive behavior. The 



problem is purely interpretation-related and concerns the 
issue of a proper Markovian limit. 

In cases with large memory effects, the formal Marko- 
vian limit, corresponding to yV(z — > 0), does not give a 
reasonable description of the system dynamics, although 
it yields correct stationary quantities, like the mean cur- 
rent. Since the noise (and also higher-order cumulants) is 
a time integral of a transient quantity, namely a current- 
current correlation function, the formal Markovian limit 
of the noise in these cases is a somewhat unphysical quan- 
tity. The problems with the interpretation of a Marko- 
vian limit also influence the interpretation of the non- 
Markovian corrections (via, e. g. negative memory times). 
Bluntly, a physically meaningful result is arbitrarily split 
into two additive parts, Markovian and non-Markovian, 
that each do not necessarily have a reasonable physical 
interpretation. The full result, however, is correct and 
physically plausible. These effects are well illustrated 
and can be understood by studying exactly solvable cases 
such as the one in Ref. ^ 

In this section, we have only briefly touched upon var- 
ious open questions and subtleties associated with in- 
terpretations of non-Markovian dynamics. However, the 
exact method developed in this paper paves the way for 
future systematic studies of memory effects in connection 
with electronic noise and counting statistics. 



VII. CONCLUSIONS 

We have presented a detailed derivation of a recursive 
scheme for evaluating high-order cumulants of transport 
through Goulomb-blockade nanostructures with many 
states and non-Markovian dynamics. In order to illus- 
trate the use of our method for Markovian systems we 
considered the counting statistics of transport through a 
two-level quantum dot and a vibrating molecule. In both 
cases, we have shown how the behavior of high-order cu- 
mulants is determined by dominating singularities of the 
cumulant generating functions. Oscillations of the high- 
order cumulants as function of the cumulant order can 
be used to locate the positions of singularities as we have 
demonstrated. We have also calculated the distribution 
of measurable currents, the so-called large deviation func- 
tion, and shown how the tails of the distributions reflect 
the high-order cumulants. In order to illustrate the use of 
our method for a non-Markovian system, we considered 
transport through a dissipative double quantum dot. For 
this system, we have studied how bath-induced dephas- 
ing affects the first three cumulants and found that effects 
of the heat bath cannot be accounted for by an effective 
detector model, when the coupling becomes strong. Fi- 
nally, we have discussed the nature and significance of 
non-Markovian dynamics in relation to counting statis- 
tics. 

The research presented in this work points to several 
interesting directions to follow. While we have focused on 
the zero-frequency current cumulants of non-Markovian 
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processes, it would be interesting to see, if the methods 
presented here could be extended to finite freouencies, 
as it was recently done for Markovian processes.cll It has 
now been firmly established that high-order cumulants of 
the counting statistics generally grow factorially with the 
cumulant order and oscillate as functions of basically any 
system parameters as well as of the cumulant order. It 
would be interesting to study in further detail how micro- 
scopic details of a system are reflected, for example, in the 
frequency of these oscillations. Such a study would shed 
new light on the information contained in high-order cu- 
mulants. Finally, we believe that the methods presented 
here will pave the way for future systematic studies of 
counting statistics in connection with non-Markovian dy- 
namics. 
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Appendix A: QR decomposition 

In this appendix we present technical details of one 
poss ible m ethod for evaluating Eqs. ( ^^ and in 
Sec. IIIB, The method we use is a standard technique 
in numerical linear algebra known as QR decomposition. 
Routines performing QR decompositions are a part of 
most common linear algebra packages such as LAPACK. 
Below we provide a piece of code implemented in MATLAB 
and subsequently explain in detail each step of the code, 
such that it can be reproduced in other programming 
languages (MATLAB itself uses an implementation of the 
QR decomposition from the underlying LAPACK library). 
The code is written in a general way but there are steps 
which are specific to the particular model considered here 
— those program lines are explicitly denoted. As a model 
system we use a double dot with 5 retained elements of 
the density matrix p and corresponding vector represen- 
tation 



with 



[1,1,1,0,0]' 



(A3) 



|p» = [Poo, Pll, Prr, Plr, Prl] (A1) 
. The trace of the density matrix can then be written 

Tr{p} = ((0|p)) = Poo + PLL + PRR (A2) 



The excerpt of the code for the evaluation of the sta- 
tionary state Eq. ( |48| ) and the pseudo-inverse (modifica- 
tion of Eq. (Il) 



(A4) 



reads: 



7o Size of the Liouville space 
N = length(W(: ,1)) ; 

7. Left zero eigenvector (MODEL DEPENDENT) 
7, Here for a double dot 
trace= [1,1, 1,0,0] ; 

7 QR decomposition with sorting the diagonal 
7 elements of 'r' in descending order 
7 (built-in routine from LAPACK) 
[q,r ,e] = qr(W) ; 

7, Consistency check - when the matrix 'W is 
7 singular the last row of the matrix 'r' 
7 should be zero 

tol = le-10; 7o Setting the tolerance 
if max(r (end, : ) ) >tol 

warning ( 'Last row of r is non-zero') 

display(r (end, : )) 

end 

7 Stationary state 

Stat = e*[r(l:N-l,l:N-l) \ r(l:N-l,end); -1]; 

7 Normalization of the stationary state 
Stat = Stat/ (trace*stat) ; 

7 Projectors 

P = kron(stat , trace) ; Q = eye(N) - P; 

7, Solution of the pseudo-inverse equation 
temp = q \ Q ; 

7 Consistency check - if the matrix 'Q' is 
7 a projector onto the regular space of 'W , 
7 the solvability condition requires the 
7 last column of the matrix 'temp' to be zero 
if max(temp(end, : ) )>tol 

warning ( 'Last row is non-zero') 

display (max (temp (end, : ) ) ) 

end 

7 Finding a particular solution of the equation 
X = e*[r(l:N-l,l:N-l)\temp(l:N-l, :) ;zeros(l,N)] ; 

7 Final fixing of the pseudo-inverse by 
7 multiplication by the projector 'Q' 
R = Q*X; 
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The code can be understood by first analyzing the 
structure of the output from the QR decomposition rou- 
tine. The QR routine takes as input a matrix W of 
dimension N x N. The output consists of three matrices 
of same dimension q, e, r, such that q • r = W • e. Here, 
q is a unitary (and thus regular) matrix, e a permuta- 
tion matrix (thus also regular) , and r an upper triangular 
matrix with decreasing diagonal elements. The column 
permutation matrix e is chosen such that abs(diag(r)) 
is decreasing. For a singular matrix W representing W, 
this implies that the last diagonal entry of r is zero and, 
therefore, the last row of r is zero. More explicitly, we 
have 



ri^2 ri^3 ri,4 
r2,2 r2,3 r2,A 
r3,4 




V 



[OllxAf-l 




ri,N \ 
r2.N 

/ 



(A5a) 



(A5b) 



with r being an upper triangular matrix with non-zero 
diagonal and dimension {N— 1) x (A^— 1) and r' a column 
vector of length A^ — 1. The QR decomposition of W 
implies for the solution |0)) to the matrix implementation 
of Eq. (E|) that 



W • = (q • r 



= 







0. (A6) 



named stat in 



Here, is a vector representation of 
the code. 

Th e block structure of the matrix r depicted in Eq. 
( A5b ) shows that (e"-'^ • 0) = c[r~^ • r', —1]^ for any num- 
ber c is a solution to the matrix equation above. We then 
findO = ce-p-i-r',-!]"^ wither l/(6-e- p-^ -r', -1]^), 
ensuring the proper normalization 0-0 = 1. We note that 
the only model dependent part of the code is the defini- 
tion of the left zero eigenvector named trace in the 
code. 

Next, we determine the pseudo-inverse R. To this end, 
we form the projectors P = ® and Q = 1 — P. For 
the Kronecker tensor produc t, th e code uses the built-in 
function kron(,). Equation (A4) can now be expressed 
in matrix form as 



W R 



(q r e 



R= Q 



R 



Q 



Since the right hand side Q of the original equation ( A4 ) 
lies in the range of W the resulting matrix equation for 
(e^^ • R) above must have a solution. This requires that 
the last row of q""'^ • Q is zero as f ollow s again from the 
block structure of r shown in Eq. ( A5b| ) (this condition 
is explicitly checked for in the code). A particular solu- 

~^ ■ [^]n-ixn 

[0]ixN 



tion of the equation is then 



with the 



rectangular matrix A being a restriction of the product 
q~^ • Q to the (A'' — 1) first rows. The pseudo- inverse is 
now fixed by multiplying this particular solution by the 
permutation matrix e and finally by the projector Q as 
follows from the discussion below Eq. (|5^). 

We have used the code on a standard PC for all of the 
examples shown in this paper and it is efficient both in 
terms of memory and CPU time for N x N matrices with 
A^ being up to several thousands. For larger matrices, the 
direct evaluation of the pseudo-inverse becomes prohib- 
ited both in terms of memory as well as CPU time and 
one should use other methods such as, e.g., the iterative 
Arnoldi scheme described in detail in App. A of Rcf. |30|. 



Appendix B: Vibrating molecule 

In this appendix we describe the model of a vibrat- 
ing molecule considered in Sec. VB and derive the cor- 
responding Markovian GME. Here we follow to a large 
extent the description of the model given in Refs. |4^, 
p^ . The molecule is operated in the Coulomb blockade 
regime, where only two charge states (m = or m = 1 
additional electrons on the molecule) participate in the 
transport. We consider spinless electrons with charge — e, 
although it would be easy to include the spin degree of 
freedom. The Hamiltonian of the molecule is 



P' 



1 



+ -mownx^ + (e - eEx)d)d, (Bl) 
2mo 2 



where E is the electric field at the position of the molecule 
with mass toq and natural oscillator frequency wq. The 
electric field is determined by the bias across the molecule 
and bias-independent contributions, e.g., image-charge 
effects. The charging energy difference between and 1 
additional electron on the molecule is denoted as e. The 
molecule is tunnel-coupled to left and right electrodes 
consisting of non-interacting fermions 



Ha — ^ £fco 



a = L, R, 



(B2) 



kept at chemical potentials fia, a — L,R, and tempera- 
ture T. Tunneling processes are accounted for by a stan- 
dard tunnel-Hamiltonian 



Ht = 



ct d-l-h.c. 



(B3) 



ka,,a=L,R 



For simplicity, we neglect any position dependence of the 
tunneling amplitudes t^^ , but it would be straightforward 
to include.E3 Finally, damping of the mechanical oscilla- 
tions are described by coupling to a bath of oscillators, 
such that the full Hamiltonian reads 



H = Hs + Ht + Hl + Hr + xVb + Hb, (B4) 



where 



(B5) 
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and 



Hb = hjjjjOjaj. 



(B6) 



The coupling to the j'th oscillator, with frequency ujj and 
corresponding creation and annihilation operators a| and 
ttj, respectively, is denoted gj. 

We treat both the coupling to the electronic leads and 
the heat bath in the weak coupling approximation and it 
thus suffices to consider the time evolution of the diago- 
nal matrix elements of the reduced density matrix of the 
charge and oscillator states of the molecule. These diag- 
onal elements correspond to the energy eigenstates of the 
isolated molecule described by Hs- The eigenstates with 
m = additional electrons on the molecule are 



\m = 0,1) = [empty) (g) \l) 
with corresponding eigenenergies 

Eqi = huJn ( I 



(B7) 



(B8) 



Here, jempty) denotes the empty charge state, while 
1^) = {a^y\0)/\/U. is the Z'th oscillator state centered at 
X = 0. The operator a'^-' lowers (raises) the oscillator 
number by 1 and |0) is the oscillator ground state. With 
m = 1 additional electron on the molecule the equilib- 
rium position of the oscillator is shifted by the distance 
d = eE /rriQUJ^. The eigenstates for the occupied molecule 
are thus 



|to = 1,0 = loccupied) (8)e'^('^*-")|0, 



(B9) 



where we have introduced the dimensionless electron- 
phonon coupling 



7 



eExo 



(BIO) 



with xq = y^h/2mou>Q. The corresponding eigenenergies 
are 

Eu^e + hajo(^l + ^^ -l^hojo- (Bll) 

In the following we denote the diagonal elements of the 
reduced density matrix by Pm,i(,n, t), where n is the num- 
ber of electrons collected in the right electrode during the 



time span [0,t]. Bath- mediated transitions between dit, 
ferent vibrational states are given by the thermal ratesc3 



Wi 



(B12) 



where K characterizes the vibrational damping rate and 
/? = l/ksT is the inverse temperature. The charge trans- 
fer rates, obtained using Fermi's Golden rule, are 



= r^,)\F,i\-fiEl;}) 



(B13) 



where / is the Fermi function, r(_|_i/_i) = r^/^j are the 
bare tunneling rates, which are assumed to be energy 
independent, i. e. 



£fej, a = L,i?. (B14) 



Moreover, we have defined 



Tpi^) TP TP 

seV 



01 



seV 



(B15) 



with Y being the symmetrically applied bias, such that 
/ii — eV 12 and /ij^ = —eV 12, and the index s indicating 
whether an electron tunneled from/to the left (s = —1) 
or right (s = -|-1) lead. Finally, the matrix elements 



Fu 



{l\e 



7(at-a)|7A 



(B16) 



are the Franck-Condon overlaps between harmonic oscil- 
lator states that have been shifted spatially with respect 
to each other due to different charge occupations, cf. Eqs. 
(B7) and (B£). In the following, we assume that the 
bias dependence of the electron-phonon coupling takes 
the form 



7 = ci + 



hujQ 



C2 



(B17) 



with ci and C2 being constants.B 

Having identified all relevant transition rates the 
Markovian master equation for the diagonal elements of 
the reduced density matrix Pm,i(n,t) reads 
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— Pm.l{^,t) 



■m' .1' -f^m^l 



l'=l±l 



s=±l,m'=0,l 



l'=l±l 



(B18) 



+ E [rt'li-^^i'Pi-m,;' t) + rSli-™,;. {(l - rn)p,,, {n-l,t) + mpo,v {n + l, t)} 



We introduce the counting fi eld v ia the transformation Pm,i{Xi't) = X^n Pni-K^i The corresponding master 
equation, obtained from Eq. ( B18D , reads 



dt 



E 



V=l±l 



l'=0 
s=±l,m'=0,l 



l'=l±l 



(B19) 



OO 

+ E +rit'li-„M' {e^Ml - m)pu'{x,t)+e-^^mpoM{x,t)} 



/'=0 



The elements Pm,iiXjt) ^-re collected in the vector p{x,t), whose equation of motion reads 



(B20) 



The matrix elements of W{x) are identified from Eq. (B19) 



Appendix C: Double dot system 



In this appendix we derive the expression for the memory kerriel-siven in Eq. (p7|). We take as our starting point 
the full Hamiltonian in Eq. (86). Following Gurvitz and PrageJi^S we project out the electronic leads in order to 
obtain an equation of motion for the reduced density matrix a = [c'm^ ^ ll^ <^ rRt ^ lR: <^ rl)^ of the three electronic 
states |0), and the bath of oscillators. The off-diagonal elements ctoq and CTao, a — L,R, between states with 

different charge occupation num bers are decoupled from the rest and can therefore be disregarded. Following the 
procedure described in Ref. 136 we find 








Tr 





\ 


Tl 








iTc 


-iTc 








-Tr 


-iTc 


iTc 





iTc 


-iTc 


-is - VrI2 








-iTc 


iTc 





le - Tr/2) 



a(t) - i 



f [HB,&Q0{t)]\ 
[HB,^LL{t) 
[HB,^RR{t) 

[HB,a-LRit) 

\[HB,ijRL{t)]J 



( \ 

~'[VB,crRR{t) 

{VB,CTLR{t)} 



(CI) 



where curly brackets denote anti-commutators {A^B} = AB + BA and we have taken h — 1. The equation is 
valid when a large-bjas is driving electrons through the double dot from the left lead to the right lead with energy- 
independent rates 



L36 



(C2) 



At this stage, the expression is valid to all orders in the tunnel coupling Tc- Due to the large-bias assumption, the 
energy Eq of the 'empty' state |0) drops out of the problem. 

We now define the electronic occupation probabilities pi = TvB{(^ii\, i = 0,L,R, where Tr^ is a trace over the 
bosonic degrees of freedom. For these probabilities we readily find 



d_ 

di 
d 

dt 
d_ 
di 



Po{t) = -TLPoit) + ^RPRit), 

Pdt) = TLPo{t) - 2TJm [TvB{(TLR{t)}] , 
pR{t) = -TRPR{t) + 2TJm [TrB{<7Li?(i)}] 



(C3) 
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We proceed by considering the equation of motion for (Jlr obtained from Eq. (CI) 



^LR{t) = -(ie + TR/2)aLR{t) - i[h'^b^ aLR{t) - aLR{t)H'g >] + iT,[aLL{t) - ^RR{t)], 



dt 



(C4) 



having defined i^jj^"* = Hb ± Vr- Its solution formaUy reads 



aLR{t) = zT, / df'e-(-+rH/2)(*-Oe-ifr (*-*') [^^^t') _ aRRit')] e'^B'it-n + e-^^'+^'^^'^^e-'"- 'aLR{0)e'"-'' 



(C5) 

The first term enters the memory kernel below, while the second term enters the inhomogeneity. In order to obtain 
a closed system of equations for the three probabilities in Eq. ( |C3| ) , we assume that the bath of oscillators between 
each tunneling event reaches a local equilibrium corresponding to the given charge state. This corresponds to the 
decoupling 



in Eq. (C5), where 



aRRit) PR{t) (E, a^-\(3) 



(C6) 



(C7) 



and (3 — l/ksT is the inverse temperature. The approximation is valid to lowest order in T"^. Note that no Markov 
approximation is made in this step. We then find 



d 
dt 



poit) = -TLPo{t) + TB,PR{t), 



dt 
d 

dt' 



ItPlW = TlPoW - / dt' \T^+\t - t')pL{t') - T^s \-t - t')PR{t')\ " lit), 



(C8) 



PR{t) = -TRPR{t)+ / dt'\T'^ + \t-t')pL{t')-T'^B\t-t')pR{t')]+^{t), 



where the inhomogeneity is of the form 7 = (0, —7, 7)"^ and the bath-assisted hopping rates are defined and evaluated 
below. An explicit expression for the inhomogeneity will not be given in this work as we are only considering the 
long-time limit for which the inhomogeneity is irrelevant. By switching to Laplace space, the memory kernel given in 
Eq. ( |87| ) is identified from Eq. ( ^8| ) after the counting field has been incorporated via the substitution — >• Trc^^ 
in the first line of Eq. (C8). In this example, the counting field enters the Markovian part of the kernel, and we do 
not encounter any of the problems described in Sec. VI B| . 



r 



In the equations above we have defined the bath-assisted 
hopping rates 



r^'it) = 2r2Re 



-(ie+rH/2)t (±) 



9'^' it) 



in terms of the bath correlation functions 



(+). 



(C9) 



(CIO) 



(±) 



These bath correlation functions cau-be evaluated using 
standard many-particle techniques. t23 First we introduce 
a polaron transformation of the form 



(Cll) 



which removes Vr from the bath correlation functions, 
since 



(C12) 



By insertion of the identity = 5^=^) [S'(=^)]t = 1 in 
Eq. ( |C10| ), we get 

g(^){t) = TrB|l(+)e-'^«^'*l(+)l(±)a(±)(/3)l(±) 
xl(-)e^<'*l(-)} 

from which standard algebra leads to 

g^^\t) = Trs{e^'*^(±*)<7(/3)e±^'^(")} 
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Here the thermal density matrix of the bath is 



/TrB{ 



(C15) 



and A{t) — e*^^*Ae '-f^^t Since Hb corresponds to free 
boson s we can write the bath correlation function of Eq. 
( |CT1 as 



with the bosonic correlation function 



W{t) 



(i^(0))o-(i(t)i(0))o 



(C16) 



(C17) 



The correlation function can be evaluated using stan- 
dard niethoijls from the field of quantum dissipatiire. 



systemsE3li2j or using a Green's functions approach.! 
Here, we just quote the final result 



W{t) 
with 



doj — ^{[1— cos(a;i)] coth 



I3uj 



-i sin(aji)}, 
(C18) 

(C19) 



being the spectral function of the heat bath. In this work 
we consider for simplicity Ohmic dissipation character- 
ized by a coupling strength a and a frequency cut-off Wc, 
such that the spectral function reads 



Jn(w) = 2acje~'^/'' 



(C20) 



For Ohmic dissipation, the correlation function is well 
known and readsllH3 



Wit) = -2a In 



|rg(l + 7? + »t//?)|^ 



(C21) 



where Te{x) is the Euler Gamma function and — 

We consider energy scales and temperatures lower than 
the cut-off Wc, such that-jf.^ 1. In that limit, Eq. ( C21 ) 
can be approximated aslH^ 



W{t) = 2a In 



sinh[z7r77(l 



ct)] 



sinh[i7r77] 



(C22) 



using only elementary functions. We can then calculate 
analytically the Laplace transform of the bath correlation 
function. Using the integral identity 



/•oo 

/ dte"^* sinh(t + x) 
Jo 



2v g-^a 
z + y 



2^1 



y + z 



(C23) 



where 2^'i(a, &, c, z) is the Gauss Hypergeometric func- 
tion, we obtain 



P\ [1 



Try 2a + Pz/iT 

f3z 



X 2F1 



a + ■ — , 2a, 1 + a 
27r 



(C24) 



vahd for a > and l//3wc, \ z\/uJc<^ 1- Without couphng 
t o th e heat bath, we would have W{t) = 0, and from Eq. 
( |C16 ) we would obtain g^^^z) — 1/z. The property of 
the Hypergeometric function, 2Fi {a, 0, c, z) = 1 for any 
value of a, c and z, shows that Eq. (C24) indeed simplifies 
to this result in the limit a — >■ 0. 

By Laplac e t ransforming the bath-assisted hopping 
rates in Eq. (|C9|), we finally find for real z 



-T![g^^Hz+) + g^^Hz^)] 



die-^*Re 



-{ie+rR/2)t(±) 



(C25) 

with z± ~ z±ie + Tii/2. Here, we have used the relation 



(026) 



for any complex z, which follows directly from the sym- 
metry property 



[W{t)]* = W{~t) 



(027) 



of the bosonic correlation function in Eq. (C18) in com- 
bination with the expression in Eq. ( |016| ). 

Before concluding this appendix, we consider the 
regime, where the coupling to the heat bath is weak and 
the bath temperature is high. Below, we show that the 
dynamics of the DQD in that regime can be described 
by a chapfte- 



rate TdB. 



Kctor model with an effective dephasing 
We derive the dephasing rate starting 
from Eq. (|C5). For weak couplings, the bath remains 
unaffected by the electronic state-of the DQD, and we 
can perform a decoupling readingl£3 



(yii{t) ~ Pi{t) ® a p. 



(028) 



Using this decoupling in Eq. (05) and tracing out the 



bath degrees of freedom we obtain 
Trs {cFLRit)} = 

iT, f dt'e-(-+i^-/2)(*-*')5(t - t') [pL{t') - PR{t')] 
Jo 

(029) 



having omitted the inhomogeneity entering Eq. (|0q), 
since we are only interested in long-time properties. The 
(single) bath correlation function g{t) is now 



5(t)=TrB{e-^^"'(*)a,e'^^"'W} 



(O30) 
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Using the polaron transformation in Eq. ( pll ), we read- 
ily find 

g{t) = (e'M0)^-■2^A{t)^^A^^)\^^ ^ g-Rc[W(t)] (^31) 

where W{t) is given in Eq. ( |Cf8D . 

For the Ohmk^bath, described by Eq. (C20), it is easy 
to demonstrateE3 in the long-time limit jit ^ h/2 that 
Re[W^(t)] « Tdt, having defined the rate 



(C32) 



Til — ^ankBT. 



Following similar steps as those leading to Eq. ( p9[ ) we 
find a bath-assisted hopping rate reading 



(C33) 



or in Laplace space 



Tb{z)^2Ti 



z + TR/2 + Td 
^e^ + {z + Tr/2 + Tdf 



(C34) 



Without coupling to the heat bath, we have F^ = 0, and 
the only broadening mechanism is the escape rate F^ of 
electrons to the right lead, which gives the hopping rate 
a width of Tr/2. With weak coupling to the heat bath 
the rate is additionally broadened by F^, and the total 
dephasing rate is Tii/2 + Td- This is similar to the results 
charge detector model with dephasing 



am 



obtained, 
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